In [1]:
# Setting up a model and a mesh for the MT forward problem

In [62]:
import SimPEG as simpeg, sys
import numpy as np
from SimPEG import NSEM
sys.path.append('/Volumes/MacintoshHD/Users/thibautastic/Dropbox/PhD_UBC/telluricpy')
import telluricpy
import matplotlib.pyplot as plt
%matplotlib inline
import copy

In [3]:
# Define the area of interest
bw, be = 557100, 557580
bs, bn = 7133340, 7133960
bb, bt = 0,480

Build the mesh

Design the tensors

hSize,vSize = 25., 10 nrCcore = [15, 8, 6, 5, 4, 2, 2, 2, 2] hPad = simpeg.Utils.meshTensor([(hSize,10,1.5)]) hx = np.concatenate((hPad[::-1],np.ones(((be-bw)/hSize,))hSize,hPad)) hy = np.concatenate((hPad[::-1],np.ones(((bn-bs)/hSize,))hSize,hPad)) airPad = simpeg.Utils.meshTensor([(vSize,13,1.5)]) vCore = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcore,(simpeg.Utils.meshTensor([(vSize,1),(vSize,8,1.3)])))])[::-1] botPad = simpeg.Utils.meshTensor([(vCore[0],8,-1.5)]) hz = np.concatenate((botPad,vCore,airPad))

Calculate the x0 point

x0 = np.array([bw-np.sum(hPad),bs-np.sum(hPad),bt-np.sum(vCore)-np.sum(botPad)])

Make the mesh

meshFor = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)


In [76]:
# Build the Inversion mesh
# Design the tensors
hSizeI,vSizeI =  25., 10.
nrCcoreI = [12, 6, 6, 6, 5, 4, 3, 2, 1]
hPadI = simpeg.Utils.meshTensor([(hSizeI,9,1.5)])
hxI = np.concatenate((hPadI[::-1],np.ones(((be-bw)/hSizeI,))*hSizeI,hPadI))
hyI = np.concatenate((hPadI[::-1],np.ones(((bn-bs)/hSizeI,))*hSizeI,hPadI))
airPadI = simpeg.Utils.meshTensor([(vSizeI,12,1.5)])
vCoreI = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcoreI,(simpeg.Utils.meshTensor([(vSizeI,1),(vSizeI,8,1.3)])))])[::-1]
botPadI = simpeg.Utils.meshTensor([(vCoreI[0],7,-1.5)])
hzI = np.concatenate((botPadI,vCoreI,airPadI))
# Calculate the x0 point
x0I = np.array([bw-np.sum(hPadI),bs-np.sum(hPadI),bt-np.sum(vCoreI)-np.sum(botPadI)])
# Make the mesh
meshInv = simpeg.Mesh.TensorMesh([hxI,hyI,hzI],x0I)
meshFor = copy.deepcopy(meshInv)

In [77]:
print np.sum(vCore)
print meshFor.nC
print meshFor


1038.9319582
99456
  ---- 3-D TensorMesh ----  
   x0: 554291.75
   y0: 7130531.75
   z0: -4530.95
  nCx: 37
  nCy: 42
  nCz: 64
   hx: 961.08, 640.72, 427.15, 284.77, 189.84, 126.56, 84.38, 56.25, 37.50, 19*25.00, 37.50, 56.25, 84.38, 126.56, 189.84, 284.77, 427.15, 640.72, 961.08
   hy: 961.08, 640.72, 427.15, 284.77, 189.84, 126.56, 84.38, 56.25, 37.50, 24*25.00, 37.50, 56.25, 84.38, 126.56, 189.84, 284.77, 427.15, 640.72, 961.08
   hz: 1393.75, 929.17, 619.45, 412.96, 275.31, 183.54, 122.36, 81.57, 2*62.75, 3*48.27, 4*37.13, 5*28.56, 6*21.97, 6*16.90, 6*13.00, 12*10.00, 15.00, 22.50, 33.75, 50.62, 75.94, 113.91, 170.86, 256.29, 384.43, 576.65, 864.98, 1297.46

In [78]:
# Save the mesh
meshFor.writeVTK('nsmesh_CoarseHKPK1.vtr',{'id':np.arange(meshFor.nC)})
nsvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh_CoarseHKPK1.vtr')

In [79]:
nsvtr


Out[79]:
(vtkRectilinearGrid)0x1192d9a70

In [80]:
topoSurf = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile('../../Geological_model/CDED_Lake_Coarse.vtp'))
activeMod = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsvtr,topoSurf)

In [81]:
topoSurf


Out[81]:
(vtkPolyData)0x1192d9ad0

In [82]:
#telluricpy.vtkTools.io.writeVTUFile('activeModel.vtu',activeMod)

In [83]:
# Get active indieces 
activeInd = telluricpy.vtkTools.dataset.getDataArray(activeMod,'id')

In [84]:
# Make the conductivity dictionary
# Note: using the background value for the till, since the extraction gets the ind's below the till surface
geoStructFileDict = {'Till':1e-4,
 'PK1':5e-2,
 'HK1':1e-3,
  'VK':5e-3}

In [85]:
# Loop through
extP = '../../Geological_model/'
geoStructIndDict = {}
for key, val in geoStructFileDict.iteritems():
    geoPoly = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile(extP+key+'.vtp'))
    modStruct = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(activeMod,geoPoly,extBoundaryCells=True,extInside=True,extractBounds=True)
    geoStructIndDict[key] = telluricpy.vtkTools.dataset.getDataArray(modStruct,'id')

In [86]:
# Make the physical prop
sigma = np.ones(meshFor.nC)*1e-8
sigma[activeInd] = 1e-3 # 1e-4 is the background and 1e-3 is the till value
# Add the structure
for key in ['Till','PK1','HK1','VK']:
    sigma[geoStructIndDict[key]] = geoStructFileDict[key]

In [87]:
fig = plt.figure(figsize=(10,10))
ax = plt.subplot(111)
model = sigma.reshape(meshFor.vnC,order='F')
a = ax.pcolormesh(meshFor.gridCC[:,0].reshape(meshFor.vnC,order='F')[:,20,:],meshFor.gridCC[:,2].reshape(meshFor.vnC,order='F')[:,20,:],np.log10(model[:,20,:]),edgecolor='k')

ax.set_xlim([bw,be])
ax.set_ylim([-0,bt])
ax.grid(which="major")
plt.colorbar(a)
ax.set_aspect("equal")



In [88]:
# Save the model
meshFor.writeVTK('nsmesh_CoarseHKPK1_NoExtension.vtr',{'S/m':sigma})

In [89]:
import numpy as np

In [90]:
# Set up the forward modeling
freq = np.logspace(5,3,16)
np.save('MTfrequencies',freq)

In [91]:
freq


Out[91]:
array([ 100000.        ,   73564.22544596,   54116.95265465,
         39810.71705535,   29286.44564625,   21544.34690032,
         15848.93192461,   11659.1440118 ,    8576.95898591,
          6309.5734448 ,    4641.58883361,    3414.54887383,
          2511.88643151,    1847.84979742,    1359.35639088,    1000.        ])

In [ ]:


In [92]:
# Find the locations on the surface of the model.
# Get the outer shell of the model
import vtk
actModVTP = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.extraction.geometryFilt(activeMod))
polyBox = vtk.vtkCubeSource()
polyBox.SetBounds(bw-5.,be+5,bs-5.,bn+5.,bb-5.,bt+5)
polyBox.Update()
# Exract the topo of the model
modTopoVTP = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(actModVTP,telluricpy.vtkTools.polydata.normFilter(polyBox.GetOutput()),extractBounds=False)
telluricpy.vtkTools.io.writeVTPFile('topoSurf.vtp',actModVTP)

In [93]:
# Make the rxLocations file
x,y = np.meshgrid(np.arange(bw+12.5,be,25),np.arange(bs+12.5,bn,25))
xy = np.hstack((x.reshape(-1,1),y.reshape(-1,1)))
# Find the location array
locArr = telluricpy.modelTools.surfaceIntersect.findZofXYOnPolydata(xy,actModVTP) #modTopoVTP)
np.save('MTlocations',locArr)

In [94]:
telluricpy.vtkTools.io.writeVTPFile('MTloc.vtp',telluricpy.dataFiles.XYZtools.makeCylinderPtsVTP(locArr,5,10,10))

In [95]:
# Running the forward modelling on the Cluster.
# Define the forward run in findDiam_MTforward.py

In [57]:
%matplotlib qt
#sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
import interactivePlotFunctions as iPf


---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-57-dab448487c8f> in <module>()
      1 get_ipython().magic(u'matplotlib qt')
      2 #sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
----> 3 import interactivePlotFunctions as iPf

ImportError: No module named interactivePlotFunctions

In [60]:
# Load the data
mtData = np.load('MTdataStArr_nsmesh_HKPK1Coarse_noExtension.npy')
mtData


Out[60]:
array([ (100000.0, 557112.5, 7133352.5, 430.0, (-1.0542164193648105-1.078716532478167j), (-23.140299596973666-19.826260505391872j), (23.389975746284172+16.906399777383363j), (0.8754499581492817+0.9556685008070969j), (-0.014494125205579434-0.07954239614803678j), (-0.041490314794175386-0.03416351101035116j)),
       (100000.0, 557137.5, 7133352.5, 430.0, (-0.8021625305031573-1.0706053468599814j), (-22.468510894286695-20.396746977494697j), (21.559482953002853+16.879166721314032j), (0.7266645327124629+1.013308657415905j), (-0.026078011261232766-0.06460099808114624j), (-0.05167463997663666-0.045362677468930454j)),
       (100000.0, 557162.5, 7133352.5, 430.0, (-1.2769734094735785-1.5464705772315357j), (-21.490658131993406-20.50650121756002j), (19.940463496713754+16.6686274356102j), (1.4318914664525015+1.704947597305036j), (-0.03226267161532167-0.0512420800212148j), (-0.06230127502660536-0.058274513337287906j)),
       ...,
       (1000.0, 557512.5, 7133952.5, 440.0, (1.197916611440613+0.6451022424541802j), (-4.228960881489435-3.790603760553026j), (7.223814873545087+5.52962654317978j), (-1.424563427490714-0.7286143613991704j), (-0.053654389267759194-0.03678133951175098j), (-0.04350593252014436-0.02697896861303861j)),
       (1000.0, 557537.5, 7133952.5, 440.0, (1.3258772950342594+0.7024381786126792j), (-4.962594275818345-4.431871833081414j), (7.141032075670906+5.614100044153096j), (-1.3994160679868353-0.7107553643963683j), (-0.06562500121402796-0.0465067927411497j), (-0.037998144645303426-0.023423112855156056j)),
       (1000.0, 557562.5, 7133952.5, 440.0, (2.0065745357923657+1.0828162579914034j), (-7.756236722381384-6.919293860657188j), (7.06987531573681+5.7449365897792015j), (-1.3988927543895338-0.7366594628536881j), (-0.06568332202732031-0.04746883772119973j), (-0.035034415429755314-0.021370552248376425j))], 
      dtype=[('freq', '<f8'), ('x', '<f8'), ('y', '<f8'), ('z', '<f8'), ('zxx', '<c16'), ('zxy', '<c16'), ('zyx', '<c16'), ('zyy', '<c16'), ('tzx', '<c16'), ('tzy', '<c16')])

In [61]:



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-61-b71e891fa06a> in <module>()
----> 1 rxT

NameError: name 'rxT' is not defined

In [110]:
iPf.MTinteractiveMap([mtData])


Out[110]:
<interactivePlotFunctions.MTinteractiveMap at 0x7f3afe78ee10>

In [104]:
# Looking at the data shows that data below 100Hz is affected by the boundary conditions, 
# which makes sense for very conductive conditions as we have.
# Invert data in the 1e5-1e2 range.

Run the inversion on the cluster using the inv3d/run1/findDiam_inversion.py


In [106]:
drecAll = np.load('MTdataStArr_nsmesh_0.npy')

In [109]:
np.unique(drecAll['freq'])[10::]


Out[109]:
array([    100.        ,     158.48931925,     251.18864315,
           398.10717055,     630.95734448,    1000.        ,
          1584.89319246,    2511.88643151,    3981.07170553,
          6309.5734448 ,   10000.        ,   15848.93192461,
         25118.8643151 ,   39810.71705535,   63095.73444802,  100000.        ])

In [93]:
# Build the Inversion mesh
# Design the tensors
hSizeI,vSizeI =  25., 10.
nrCcoreI = [12, 6, 6, 6, 5, 4, 3, 2, 1]
hPadI = simpeg.Utils.meshTensor([(hSizeI,9,1.5)])
hxI = np.concatenate((hPadI[::-1],np.ones(((be-bw)/hSizeI,))*hSizeI,hPadI))
hyI = np.concatenate((hPadI[::-1],np.ones(((bn-bs)/hSizeI,))*hSizeI,hPadI))
airPadI = simpeg.Utils.meshTensor([(vSizeI,12,1.5)])
vCoreI = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcoreI,(simpeg.Utils.meshTensor([(vSizeI,1),(vSizeI,8,1.3)])))])[::-1]
botPadI = simpeg.Utils.meshTensor([(vCoreI[0],7,-1.5)])
hzI = np.concatenate((botPadI,vCoreI,airPadI))
# Calculate the x0 point
x0I = np.array([bw-np.sum(hPadI),bs-np.sum(hPadI),bt-np.sum(vCoreI)-np.sum(botPadI)])
# Make the mesh
meshInv = simpeg.Mesh.TensorMesh([hxI,hyI,hzI],x0I)

In [94]:
meshInv.writeVTK('nsmesh_HPVK1_inv.vtr',{'id':np.arange(meshInv.nC)})

In [101]:
nsInvvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh_HPVK1_inv.vtr')
activeModInv = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsInvvtr,topoSurf,extBoundaryCells=True)

In [102]:
sigma = np.ones(meshInv.nC)*1e-8
indAct = telluricpy.vtkTools.dataset.getDataArray(activeModInv,'id')
sigma[indAct] = 1e-4

In [103]:
meshInv.writeVTK('nsmesh_HPVK1_inv.vtr',{'id':np.arange(meshInv.nC),'S/m':sigma})

In [108]:
from pymatsolver import MumpsSolver


---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-108-d6799536a06f> in <module>()
----> 1 from pymatsolver import MumpsSolver

ImportError: cannot import name MumpsSolver

In [106]:
pymatsolver.AvailableSolvers


Out[106]:
{'DiagonalSolver': True,
 'Mumps': False,
 'TriangleFortran': False,
 'TrianglePython': True}

In [116]:
NSEM.Utils.skindepth(1000,100000)


Out[116]:
50.329212104487041

In [117]:
np.unique(mtData['freq'])[5::]


Out[117]:
array([   1000.        ,    1584.89319246,    2511.88643151,
          3981.07170553,    6309.5734448 ,   10000.        ,
         15848.93192461,   25118.8643151 ,   39810.71705535,
         63095.73444802,  100000.        ])

In [ ]: